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We study numerically the interplay of phase coherence and vortex-glass state in two-dimensional 
Josephson-junction arrays with average rational values of flux quantum per plaquette / and random 
dilution of junctions. For / = 1/2, we find evidence of a phase coherence threshold value x s , below 
the percolation concentration of diluted junctions x p , where the superconducting transition vanishes. 
For x s < x < x p the array behaves as a zero-temperature vortex glass with nonzero linear resistance 
at finite temperatures. The zero-temperature critical currents are insensitive to variations in / in 
the vortex glass region while they are strongly / dependent in the phase coherent region. 
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The study of the structure of vortex-lattice states in 
disordered superconductors in a magnetic field have at- 
tracted much recent interest. In three-dimensions a true 
vortex-glass transition!!! at finite temperature is possi- 
ble for strong disorder in unscreened superconductorsH, 
while for weak disorder a Bragg-glass phase with quasi- 
long range order has been prcoiosedH. In two-dimensions, 
however, vortex-glass modelsua and experiments on su- 
perconducting filmsH show that vortex-glass order and 
phase coherence are destroyed at any finite tempera- 
ture with a nonzero but exponentially small resistivity 
in the large disorder limit whereas in the weak disor- 
der limit the situation is less clear and may depend on 
the particular model of disordered superconductor. Ran- 
dom diluted Josephson-junction arraysJaave been used 
to model disordered superconductorsQllJ and can also 
be fabricated path controlled amount of disorder in two- 
dimensionsE3~til. For a regular array in perpendicular 
magnetic field with a rational flux quanta per cell /, the 
ground sate consists of a periodic pinned vortex lattice, 
with additional discrete-symmetries resulting from com- 
mensurability effectsHE3, and phase coherence and vor- 
tex order is possible. Thus, diluted arrays in a magnetic 
field can provide a convenient experimental model system 
to investigate the effects of weak and strong disorder on 
initially pinned vortex lattices and the interplay of phase 
coherence and vortex glass states in two dimensions. In 
particular, in order to understand transport properties 
near percolation threshold in recent experimentsEJ on di- 
luted arrays in a magnetic field, it is important to know if 
disorder and temperature fluctuations can destroy phase 
coherence at long length scales and [the nature of vortex 
order in this regime. A recent studyEH of a model of ran- 
dom Josephson-junction arravs.with a particular type of 
disorder (positional disorder in a magnetic field sug- 
gests that no transition is possible even for weak disorder 
in the thermodynamic limit but it is not clear if this sce- 



nario would apply in general. In fact, random dilution 
does not explicitly introduces random phase shifts across 
the junction unlike positional disorder. In addition, an 
earlier study of the ground-state stability of a diluted ar- 
ray showali that, in presence of a magnetic field with an 
average rational value of flux quantum per plaquette /, 
phase coherence is possible at nonzero temperatures and 
the transition temperature only vanishes, for increasing 
dilution of junctions x, at a critical value x s (f) below 
the geometric percolation threshold x p where the transi- 
tion would vanish in the absence of the magnetic field. 
A vortex-glass phase with zero-temperature transition 
should also appear at a critical concentration x v > x s , 
with nonzero resistivity at finite temperatures but hav- 
ing a diverging short-range correlation length £ cx T~ v 
which is expected to determine the nonlinear behavior 
of the current voltage characteristics. An upper bound 
for the phase-coherence region is set by the behavior at 
/ = 1/2 since higher order rational values are expected to 
be much less stable with a corresponding threshold value 
which may be too small to detect numerically. These re- 
sults rely on the finite-size behavior of defect energy in 
the ground state which are inaccessible experimentally. 
However, experiments often measure transport proper- 
ties and it is of great interest to know how these effects 
could show up in the behavior of the current- voltage char- 
acteristics. 

In this work, we present the results of extensive dy- 
namical simulations of the current- voltage characteristics 
of resistively shunted Josephson-junction arrays with an 
average flux quantum per plaquette / and random dilu- 
tion of junctions. We find evidence of a phase coherence 
threshold value x s < x p as indicated in Fig. 1. For 
x < x s , the superconducting transition occurs at finite 
temperatures while for x s < x < x p the array behaves as 
a zero-temperature vortex glass with nonzero and ther- 
mally activated linear resistance at finite temperatures 
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and diverging short-range correlation length ( a T _ ". A 
current-voltage scaling analysis provides an estimate of 
v ~ 2. In the vortex-glass region, the zero-temperature 
critical currents are roughly insensitive to changes in /. 

We consider a two-dimensional array of superconduct- 
ing grains coupled to its nearest neighbors by resistively 
shunted Josephson_iunctions and with current conser- 
vation at each sitetj. The equations of motion for the 
phases 9i of the superconducting order para meter located 
at site i of the lattice can be written asEiEj 



^2(9i - 6j) = - Y^[Iij sin(6i - 6j - A l3 ) + nq] 

3 3 

(1) 



2eR D 



where R is a uniform shunt resistance, T)ij(t) is 
a thermal noise with correlations {f]ij(t)ilkl('t')) — 
2ksT ' j ' R 5ij,kiS(t — t') and Zy is the junction critical cur- 
rent. The bond variables correspond to the line 
integral of the vector potential and are constrained to 
yV ■ Aij — 2nf, about each elementary plaquette of the 
reference (undiluted) lattice. For simplicity we consider a 
square lattice array and bond dilution of junctions. The 
qualitative behavior and critical exponents presented be- 
low should remain the same for other choices of dilution 
and for triangular arrays. Dilution of junctions is intro- 
duced by taking 1^ ~ with probability x and Ty = 7 , 
a constant, with probability 1 — x. Dimcnsionless quanti- 
ties are used with time in units of r = ti/2eR J , current 
in units of I , voltages in units of R I and temperature 
in units of HI /2ekB- A total current I is imposed urik 
formly in the array using periodic boundary conditionscJ 
with current density J — I / L, where L is the system size 
and the average electric field E is obtained from the volt- 
age V across the system as E — V/L = (%/2e)(Qi — Oj). 
We use periodic boundary conditions in order to elimi- 
nate possible edge contributions to the resistance due to 
diluted junctions near the boundary which could arise 
from open boundary conditionsEU. System sizes ranging 
from L — 8 to L = 128 were used in the calculations 
with a time step At = 0.0 7r and the results averaged 
over 10 to 500 random diluted configurations of junctions 
depending on the system size. 

Fig. 1(b) shows the behavior of the critical current 
density J c where a nonzero voltage appears at zero tem- 
perature. At low values of x, the behavior of J c strongly 
depends on the rational frustration / = p/q, as indicated 
for / = 1/2 and / = 1/4 in the Figure, but becomes 
roughly insensitive to / for dilutions larger than a critical 
value much below the percolation threshold a;™—. This is 
consistent with the proposed vortex glass phaseO for the 
range x v (f) < x < x p where dilution threshold be- 

low which vortex-lattice order remains. In the undiluted 
case and for small enough x the ground state energy and 
critical currents correlate with the ordering xtf q x q unit 
cells and so are very sensitive to the q valueEJj. However, 
for x > x v , vortex-lattice order is completely destroyed 
at long-length scales and its stability and therefore the 
critical current should be less sensitive to q. Since it is ex- 



pected thatS 2^(1/4) < x v (l/2) and x v > x s , the change 
in the behavior of J c in Fig. 1(b) allows for a very rough 
estimate of the phase coherence threshold for / = 1/2 as 
the value of dilution where the two curves overlap within 
the estimated crrorbars, x s (l/2) ~ 0.20(5) . 

We turn now to thermal fluctuation effects. Fig. 2(a) 
shows the temperature dependence of the nonlinear re- 
sistivity E/J at a value of dilution x = 0.1 below the 
phase-coherence threshold x s estimated above, for the 
largest systems sizes L — 64 and L ~ 128. As can be seen 
from the Figure, the linear resistivity Rl = lim./^o E/J, 
estimated from the ratio E/J when J — > 0, tends to a fi- 
nite value at high temperatures but extrapolates to very 
low values at lower temperatures, independent of system 
size, consistent with the existence of a finite temperature 
superconducting transition in the range T c = 0.3 to 0.4. 
This is confirmed by a scaling analysis of the nonlinear re- 
sistivity according to whichta measurable quantities scale 
with the diverging correlation length £ near the transition 
temperature. If the transition occurs at a finite temper- 
ature, the relaxation time diverges as £ z , where z is the 
dynamical critical exponent, and the nonlinear resistivity 
satisfy the scaling form 



(2) 



in two dimensions, where the + and — correspond to the 
behavior above and below the transition, respectively. 
For a transition in the Kosterlitz-Thouless (KT) univer- 
sality class, the correlation length should diverge expo- 
nentially as £ oc exp(b/\T/T c — ll 1 / 2 ), while for a con- 
ventional transition a power-law behavior is expected 
£ oc \T/T C — l| _ly , with an exponent v depending on 
the discrete symmetry of the pinned vortex lattice. A 
scaling plot according to Eq. (|^) can be used to verify 
the scaling arguments and the assumption of finite tem- 
perature equilibrium transition. This is shown in Fig. 
2(b), in the temperature range closest to the apparent 
T c and smallest current densities, assuming the corre- 
lation length £ has an exponential divergence as in the 
KT universality class and using b, T c and z as adjustable 
parameters so that the best data collapse is obtained. 
As shown in the Fig. 2(b), the two largest system sizes 
L = 64 and L = 128 give the same data collapse and 
so finite size effects, ignored in the scaling form of Eq. 
(^) , are not dominant for this range of temperatures and 
current densities. We estimate a transition temperature 
T c = 0.33(2) and dynamical exponent z — 1.2(2). Al- 
though this estimate is based on a scaling analysis of 
the nonlinear current-voltage characteristics, which is a 
nonequilibrium property, we find that the finite size be- 
havior of the linear resistance at T c is consistent with 
this analysis. In a finite system the divergent correlation 
length £ is cut off by the system size L at the transi- 
tion. From Eq. (|^), the linear resistance at T c should 
then scale as Rl oc L~ z . The linear resistance can be 
obtained from the Kubo formula of equilibrium voltage 
fluctuations as R L = (1/2T) / dt < V(t)V(0) >, without 
finite current effects, and can also be determined from 
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the long-time fluctuations of the phase difference across 
the systemQ. Fig. 3, shows the finite size behavior of Rl 
at different temperatures. Above our estimate of T c it 
remains finite for increasing L whereas below T c it ap- 
pears to extrapolate to zero. Right at T c , a power-law fit 
gives z = 1.07(5) which is consistent with the estimate 
from the current-voltage scaling and suggests therefore 
that the transition corresponds to the underlying equilib- 
rium behavior, ft should be noted that for the pure KT 
transition a dynamical exponent z = 2 is expected, in- 
dependent of the particular dynamics. Indeed, for / = 
and x — the same power-law fit gives z = 2.0(1) at the 
critical temperature. However, for / = 1/2, where an ad- 
ditional Ising order parameter is present, it is found that 
even for the undiluted system z < 2 using the present 
dynamicsE3. We also note that attempting a scaling plot 
using the conventional power-law correlation length gives 
a very large value for v which suggests that the exponen- 
tial form is the appropriate one. However at a; = 0, 
both forms of correlation length gives reasonable data 
collapse as expected for / = 1/2 from the single transi- 
tion scenario where the superconducting transition and 
vortex-lattice disordering transition occurs at the., same 
temperature or else at very close temperaturesESEJcj. 
Since the undiluted array at / = 1/2 is expected to have a 
transition combining the KT and Ising universalities, our 
results suggest that for -C x < x s , the superconducting 
transition is in the KT universality class (statics) while 
the vortex-lattice disordering transition of Ising symme- 
try may occur separately in presence of weak disordered. 
However, our above scaling analysis based on the diverg- 
ing phase-coherence correlation length does not allow us a 
determination of the vortex-lattice disordering transition 
since it is expected to occur within the normal phaseEij. 
Using the above scaling analysis, the transition temper- 
atures for the different dilutions can be obtained as in 
Fig. 1(a). For values of x ^ 0.1 limited data was used 
and the results are only rough estimates of the transi- 
tion temperatures. Nevertheless, the critical temperature 
as a function of dilution reasonably extrapolates to the 
threshold x s estimated from the behavior of J c in Fig. 
1(b) as discussed above. 

In contrast, for a dilution above the phase-coherence 
threshold x > x s , the linear resistance Rl is finite for 
all temperatures in the same range as indicated in Fig. 
4a. Although we can not exclude a transition at much 
lower temperatures based on these data, the behavior is 
consistent with a superconducting transition and vortex 
order occurring only at zero temperature, .as for a vortex 
glass with a zero-temperature transitioncm. This is con- 
sistent with defect energy calculations which show that 
low-energy excitations above the ground state decreases 
with system size in this range of dilutionsllll. In fact, 
Rl decreases rapidly with decreasing temperature and 
for increasing J there is a smooth crossover to nonlinear 
behavior at a critical current J n i which also decreases 
with decreasing temperature. From the well-known scal- 
ing argumentsa leading to Eq. (Q), if the transition hap- 
pens only at zero temperature then £ oc T~ v and since 



the current density scale as J oc kT/£, the crossover to 
nonlinear behavior sets in at J n i oc T 1+v which depends 
strongly on the yet unknown critical exponent v. Also, 
the linear resistivity Rl is finite at any nonzero tempera- 
ture but thermally activated, Rl oc exp(—Eb/kT). Thus 
the relaxation time r oc 1/Rl diverges exponentially for 
decreasing temperatures. We can then consider the be- 
havior of the dimensionless ratio E / ' JRl which must sat- 
isfy the scaling forrrfl 

E/JRl = g(J/T 1+ ») (3) 

if the assumption of a zero-temperature transition is cor- 
rect. In Fig. 4(b) we show the scaling plot according to 
Eq. (|J) for the lowest temperatures and current densities 
which verifies the scaling assumption and provides an es- 
timate oiv = 2.1(2) and an energy barrier Ef, = 1.2. This 
value of v is consistent with the estimate ~ 1.9 based 
on the previouS|-ijinite-size scaling of defect energy in 
the ground stateO. Similar analysis at different dilution 
x = 0.45 gives v = 2.2(2) and Ef, = 0.9 and at different 
frustration / = 1/4 gives v = 2.3(2) and E b = 1.45. We 
note that our estimate of v ~ 2 is roughly the same as the 
value obtained for the gauge-glass model offStrongly dis- 
ordered two-dimensional superconductorsura which may 
suggest a common universality class. However, it should 
be noted that, for / = 1/2, the system has a global reflec- 
tion symmetry (9i —> —9%) in addition to the rotational 
symmetrynS and one would expect, similarly to the XY 
(chiral) spin glassEZI which shares the same feature, two 
different divergent correlation lengths £ s and £ c with cor- 
responding distinct exponents v s and v c . In fact, for 
the chiral glass model a different universality class with 
v s ~ ljJias been found from a current-voltage-scaling 
analysisa. On the other hand, an analytic studyEa of the 
XY spin glass for a particular distribution of disorder 
find a common exponent. Our estimates suggest that 
this could also be the case for the present percolative 
type of disorder or else the exponents are too close to 
be resolved within the accuracy of our estimate. An ap- 
parent common universality class of vortex glass models 
with clearly distinct symmetries has also been found in 
three dimensionsEl. In addition, close to the percolation 
threshold x p , the above scaling analysis based on a single 
diverging length scale is not valid, as one must also take 
into account the percolation correlation length £ p an d t h e 
fractal nature of the system at smaller length scalestiaEJ. 

In summary, we have studied the interplay of phase co- 
herence and vortex-glass state in two-dimensional diluted 
Josephson-junction arrays with average rational values of 
frustration. For / = 1/2, we found evidence of a phase 
coherence threshold value x s much below the geometric 
percolation threshold x p . ThisJs in contrast with the con- 
clusions of Gupta and TeiteO for a Josephson-junction 
array with positional disorder where no phase coherence 
is expected at finite temperatures even for small disor- 
der at length scales much larger than a disorder depen- 
dent length. Further work is required to verify whether 
the present study only reflects the finite-length scale of 
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the system sizes used in the calculation or is a conse- 
quence of different type of disorder. In addition, since 
/ = 1/2 has a particular reflection symmetry, which is 
preserved in presence of random dilution, the behavior 
for other values of / could he-qualitatively different. On 
the other hand, experimentaij are often done on systems 
sizes comparable to our largest system size and thus the 
current- voltage scaling behavior discussed here should be 
observable. In the range x s < x < x p the array behaves 
as a zero-temperature vortex glass with activated nonzero 
linear resistance at finite temperatures and critical cur- 
rents much less sensitive to variations in / than in the 
phase coherent region. Our results suggest that the phase 
coherence threshold can be identified experimentally as 
the change in the transport properties from the weak 
to the strong disorder regime. However, the numerical 
estimates of critical quantities from the current-voltage 
scaling analysis should be regarded as rough magnitudes 
which can be measured experimentally and used to verify 
the prediction of a change in behavior of the transport 
properties in the different phases. Equilibrium simula- 
tions are required to confirm the observed critical be- 
havior and obtain more accurate estimates of the critical 
exponents. 
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FIG. 1. (a) Phase diagram of a diluted Josephson-junction 
array as a function of temperature T and concentration x of 
diluted junctions, for an average rational frustration / = 1/2. 
The superconducting phase is denoted by S, the normal phase 
by N and the short-range vortex glass state by VG. The ge- 
ometrical percolation threshold is indicated by x p and the 
phase-coherence threshold by x s . (b) Critical current den- 
sities J c as a function of dilution x for different values of 
frustration /. Critical temperatures (open circles) in (a) were 
obtained from current-voltage scaling analysis and the phase 
coherence threshold x 3 (filled circle) was inferred from the 
change in the behavior of J c (/) in (b). 
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FIG. 2. (a) Nonlinear resistivity E/J as a function of 
temperature for a dilution x — 0.1 below the phase coherence 
threshold x s and system sizes L = 64 and L — 128 (symbols 
connected by lines), (b) Scaling plot of the data (not indicated 
in (a)) for the smallest range near T c and smallest current 
densities. Open symbols correspond to L = 64 and filled ones 
to L = 128. 
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FIG. 3. Linear resistance as a function of system size L 
for different temperatures at x — 0.1. A power-law fit at 
T c = 0.33 gives an estimate of the dynamical exponent z. 



FIG. 4. (a) Nonlinear resistivity E/J as a function of tem- 
perature for a dilution x = 0.35 above the phase coherence 
threshold x s . (b) Scaling plot of the data in (a) for the low- 
est temperatures and current densities according to a T = 
transition. 
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